Thermalization and free decay in Surface Quasi-Geostrophic flows 
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We derive statistical equilibrium solutions of the truncated inviscid surface quasi-geostrophic 
(SQG) equations, and verify the validity of these solutions at late times in numerical simulations. 
The results indicate the pseudo-enstrophy thermalizes while the pseudo-energy can condense at the 
gravest modes, in agreement with previous indications of a direct cascade of pseudo-enstrophy and 
an inverse cascade of pseudo-energy in forced-dissipative SQG systems. At early times, the trun- 
cated inviscid SQG simulations show a behavior reminiscent of forced-dissipative SQG turbulence, 
and we identify spectral scaling laws for the pseudo-energy and pseudo-enstrophy spectra. More 
importantly, a comparison between viscous and inviscid simulations allows us to identify free-decay 
laws for the pseudo-enstrophy in SQG turbulence at very large Reynolds number. 
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I. INTRODUCTION 

Classical Gibbs ensemble methods have been exten- 
sively applied to Galerkin representation of turbulent sys- 
tems, and many examples can be found in the literature. 
The first studies of the statistical mechanics of discrete 
distributions of vortices in a two-dimensional (2D) flow 
using a Hamiltonian formalism can be found in [l| . Later, 
it was shown 0] that Gibbsian statistical mechanics can 
be applied to Galerkin truncations of the hydrodynamic 
and magneto hydrodynamic (MHD) equations. This al- 
lowed many studies of continuous vorticity distribution 
in truncated 2D inviscid flows 0-Q , where absolute equi- 
librium spectra for the quadratic conserved quantities of 
the system were derived. Following 0]i in the absolute 
inviscid equilibrium ensemble for three-dimensional (3D) 
Euler flows was considered. Later, 3D inviscid magneto 
hydrodynamic MHD equilibrium solutions were investi- 
gated in [1], and recently the approach was extended to 
Hall-MHD in Q. Other problems studied in this frame- 
work include geophysical flows [l^, fast rotating flows 
[ill . [T2j . and formulations of one- and two-layer quasi- 
geostrophic models [T^. In general, classical Gibbs en- 
semble methods can be applied to systems with quadratic 
conserved quantities and that satisfy the Liouville theo- 
rem, which implies the incomprcssibility of the flow in 
phase space (e.g., the space of complex Fourier modes). 

In the absolute equilibrium, fields have Gaussian statis- 
tics and the quadratic conserved quantities of the system 
have associated temperatures that can be positive or neg- 
ative. In the former case, the quantity is said to "ther- 
malize" , and the quantity in the equilibrium is equally 
distributed among all modes in the system. In the lat- 
ter case the quantity "condenses," and is accumulated 
at the gravest modes. However, forced-dissipative tur- 
bulent systems are far from equilibrium, as the effect of 
viscosity prevents relaxation towards a true equilibrium 
state, giving rise to solutions with non-zero flux (cas- 
cades) and making direct comparison with equilibrium 
solutions inadequate. In spite of this, Gibbs ensemble 



methods proved useful to predict the direction of the 
cascades in many forced-dissipative systems, depending 
on whether the quantity of interest thermalized or con- 
densed in the associated truncated inviscid system (see, 
e-g-, |1, 0, [IS [3 : although note not all inverse cascades 
in forced-dissipative systems are associated with conden- 
sates in the statistical equilibrium of the equivalent trun- 
cated ideal system [TTI. [l2jV 

The recent finding of transients and quasi-stable states 
in truncated inviscid flows, in which non-thermalized 
modes behave as in the viscous case (see, e.g., [TB - ITsj ). 
has renewed interest in classical Gibbs ensembles, indi- 
cating more information than just the direction of the 
cascades can be extracted from these systems. In par- 
ticular, it was found that at early times and as the sys- 
tem evolves towards equilibrium, a comparison between 
the inviscid system and a viscous turbulent flow can be 
achieved by considering in the inviscid system the net 
effect of the modes that have already thermalized as an 
effective viscosity acting on the non-thermalized modes 
[TgI . ITsI . This viscous- like dynamics was reported in 
many systems, including the 3D truncated Euler equa- 
tion O, Ulllil, 3D truncated rotating flows [13], 2D 
truncated MHD [T^l, and turbulent Bose-Einstein con- 
densates using the truncated Gross-Pitaevskii equations 
[20| . In this latter case, the truncated equations allowed 
the study of thermalization of Bose-Einstein condensates 
at finite temperature. 

Classical Gibbs ensembles and quasi-stable states were 
also used to study the behavior of atmospheric models 
(see, e.g., [ITI), as the need to study atmospheric and 
oceanic flow dynamics has led to a variety of approximate 
models derived from the 3D Navier-Stokes equation for 
stratified and rotating flows (see for example (22|j23 |). 
In the often used geostrophic approximation, the vertical 
component of the velocity field is assumed zero (uz = 0), 
and hydrostatic balance is solved in that direction while a 
linear balance between the Coriolis force and the pressure 
gradient is solved on the horizontal plane. Another fam- 
ily of models is given by the so-called quasi-geostrophic 
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(QG) models, which are a first order departure from the 
linear geostrophic balance on the horizontal plane. That 
the truncated QG systems satisfy the Liouville theorem, 
thus allowing the use of Gibbs ensemble methods, can be 
seen, e.g., in Ref. [2]| . 

A particular case of the QG models is the Surface 
Quasi-Geostrophic (SQG) approximation [25l - l30| . which 
describes rotating stratified flows with constant poten- 
tial vorticity. In this model, the vertical gradient of the 
stream function ip matches a scalar field (the density 
field p or the potential temperature T) at a flat surface 
z = 0. The scalar field is identified with the horizon- 
tal Laplacian (— A)^/^';/), and the equation for the ad- 
vection of this scalar by the incompressible surface flow 
u = f X VV' = {—dyijj, dxip, 0) is solved. Contrary to 3D 
QG turbulence whose dynamics is driven by large-scale 
interior potential vorticity gradients, SQG flows are en- 
tirely driven by density or potential temperature vari- 
ations at the surface. Recently, this system has been 
used to study the dynamics of the upper troposphere 
[3ll - l33| , and of the upper oceanic layers with relative ac- 
curacy down to 500 meters [33 - f36| . Turbulence in forced- 
dissipative SQG systems (including direction of cascades 
and scaling laws) has also been studied in numerical sim- 
ulations ^5, 28, 30,]. 

Besides potential applications in ocean and atmo- 
spheric dynamics, the SQG equations have interest from 
the mathematical and physical points of view. Although 
the SQG equations are 2D in nature, their dynam- 
ics share similarities with the 3D Euler equations, and 
whether a singularity in the gradients of the scalar field 
develops at finite time from smooth initial conditions 
has been the focus of many mathematical and numerical 
studies [srl - lioj . Although there is no clear-cut answer 
for the most general case, a singularity has been ruled 
out for the case in which isolevels of the scalar contain 
a hyperbolic saddle [4l| . In the same context, numerical 
comparisons between the dynamics of the inviscid and 
viscous SQG equations were carried in [l^. More re- 
cently, cascades in the SQG equations have received at- 
tention from the physical point of view, as it was found 
that large-scales in SQG flows have conformal invariant 
properties [43j . 

In this paper we investigate the classical Gibbs en- 
semble solution of the truncated inviscid SQG equations, 
and compare these solutions and their viscous-like tran- 
sient with solutions of the dissipative SQG equations. 
In particular, we derive solutions for the pseudo-energy 
and pscudo-enstrophy equilibrium spectrum, verify nu- 
merically the convergence of the truncated SQG solu- 
tions to the statistical equilibrium solutions, and com- 
pare inviscid and viscous numerical results. The statis- 
tical equilibrium solution indicates the pseudo-enstrophy 
thermalizes while the pseudo-energy condenses in SQG, 
in agreement with previous numerical results where a 
direct cascade of pseudo-enstrophy and an inverse cas- 
cade of pseudo-energy were reported [25|, H^]. In the 
viscous-like transient, we observe the development of in- 



ertial ranges with ^ k~^/^ scaling for the spectrum of the 
pseudo-enstrophy, unlike previous inviscid simulations 
where steeper power laws were observed [42|. Finally, 
we present an analogy between the non-thermalizcd frac- 
tion of the pseudo-energy and pseudo-enstrophy in the 
inviscid truncated runs, and the free decay of the pseudo- 
energy and pseudo-enstrophy in the viscous SQG equa- 
tions at very large Reynolds number. The analogy allows 
us to identify possible asymptotic behavior for very large 
Reynolds numbers. 

The paper is organized as follows. In Sec. |TT]we intro- 
duce the SQG equations and derive the statistical equilib- 
rium solutions using the canonical distribution function 
formalism. In Sec. lIIII we present a set of numerical simu- 
lations of the truncated inviscid SQG equations. We first 
compare them with the theoretical equilibrium solutions, 
and we then introduce a viscous SQG system comparing 
the time evolution of viscous and inviscid simulations. 
Finally, we present our conclusions in Sec. IIVI 

II. SQG AND THE GIBBS ENSEMBLE 

The system considered in this paper is an incompress- 
ible 2D SQG flow. Its equations are usually presented in 
the literature as part of a family of equations governing 
the advection of a scalar [2^ 

g = (-Ar/v. (1) 

This family includes the SQG equations {a — 1, which 
is the case in the present study) [2^, H^, the vorticity 
equation in an Euler flow (a = 2), and the equation of 
motion for a shallow flow in a rotating domain driven by 
uniform internal heating (a = 3) [131 . The case a = —2 
corresponds to a shallow-water QG equation in the limit 
of large length scales compared to the deformation scale 
(see [46j). 

The flow in these models is described by a stream func- 
tion ip, and governed by the so-called a-turbulence equa- 
tions which, in their inviscid unforced version, are usually 
written in the general form 

dtq + J{i',q) = 0, (2) 

where J is the Poisson bracket 

J{A, B) = d^AdyB ~ d.BdyA. (3) 

For SQG, Eq. ([1]) in Fourier space reduces to qCk) = 
|k|V)(k), and Eq. ^ reads 

dt'4' = 1^ (d^^q dy^p - d.^'ip dyq^ , (4) 

where the hats denote Fourier transformed. Inviscid SQG 
dynamics possesses two quadratic conserved quantities, 
analogous to the energy and enstrophy in 2D Euler flows, 
which are defined respectively as 

E ^ -\ I qipdxdy, (5) 
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and 



G = ^ / dxdy, 



(6) 



with A the total area of the integration domain. We 
wih refer to these magnitudes as the pseudo-energy and 
pseudo-cnstrophy respectively. Note that sometimes 
these quantities are defined with a factor 1/2 in front, 
in analogy with the energy and enstrophy in the 2D Eu- 
ler equations (see, e.g., (l^)- For the sake of simplicity, 
in the following we use the convention in [2^, [2^ and 
drop those factors. 

Note in the forced-dissipative case, Kolmogorov- 
Batchclor-Kraichnan phenomenology predicts two iner- 
tial ranges for these quantities [1^ : an inverse cascade of 
pseudo-energy with spectra E{k) ^ and G{k) ~ 
and a direct cascade of pseudo-enstrophy with spectra 
E ^ k^^/^ and G fc"^/'^. We will consider this case in 

secnni 

Writing k ~ |k|, the pseudo-energy and pseudo- 
enstrophy of each Fourier mode k arc easily related by 



G(k) = fc£;(k) = |M(k)| 



(7) 



Since Wx(k) = — ifcj,'i/'(k), Uj,(k) — ikx^'i}^), and the re- 
lation between the pseudo-energy in each mode and the 
stream function is 



A^^{kl + klM^^k^\i;\\ 



(8) 



the pseudo-energy and pseudo-enstrophy can then be 
written in Fourier space as follows 



E= 

G= fc'i^wr 



(9) 



(10) 



where N' = {1, N /2} runs over all degrees of freedom 
of the system. Note that in Eqs. ^ and ([TU| we already 
truncated the system up to a finite number of modes. 

From these relations we can derive the statistical equi- 
librium solutions for the inviscid truncated SQG system. 
The generalized Gibbs canonical distribution is 



z 



where 



Z = 



I 

V phase spa 



-(/3E+7G) 



(11) 



(12) 



is the partition function and = J^^. dui(k)du2(k); Ui 
and U2 are defined such that u(k) = Ui(k) -|-m2(k) and 
are constrained by the incomprcssibility condition V • u = 
(or k • u = in Fourier space) by 



kyitly 



0, 



(13) 



kxU2x + kyU2y = 0. 

From Eqs. and ([TU| . 

|u(k)p 

k/fceW 



G= ^ |u(k)p, 

k/fceiV' 

and the partition function results 



(14) 



(15) 



(16) 



Z 



k/fceA" 



n 2- 

k/fceAT' 



,-|ui(k)p(f+7) 



Ui(k)|dwi(k) X 



X 27r / e-l"^('')l'(^+^)|u2(k)Mii2(k). (17) 



Note we used isotropy of the velocities Ui and U2 to 
convert the first integral in the entire phase space into two 
one-dimensional integrals with respect to the absolute 
values ui = |ui| and U2 = |u2| respectively. 
Using 



2(7fc + /3)' ^'''^ 



which holds provided i?e(7 + /3/k) > 0, the partition 
function results 



Z 



n 

k/keN' 



nk 



■yk + 13 



n 

k/kcN' 



(19) 



The most probable probability density function can be 
obtained by other methods. A Lagrange multiplier 
method can be used (see [2l|). Another elegant way is 
to use properties of n-dimensional normal distributions 
(see, e.g., Q)- All these methods yield the same result. 

From the partition function we can derive the mean 
modal intensity spectra of pseudo-energy and pseudo- 
enstrophy 



G(k) 



91nZk 

d(3 



i91nZk 
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7fc + /3 



2k 



■yk + /3 



(20) 



(21) 



These spectra give the pseudo-energy and pseudo- 
enstrophy in each mode k, and are functions of its mod- 
ulus only. The usual isotropic spectrum is obtained by 
integratiiig Eqs. ([20|) and (f2T|) . resulting, e.g., E{k) = 
7rfci?(k) [5| provided the modes are dense enough over 
the entire spectrum. We then finally get the expressions 
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for the isotropic spectra of pseudo-energy and pseudo- 
enstropliy, 

/2-wh 
E(k)kdip = TTkE{k) = (22) 



1 .000 F 



G{k) = / G(k)kdip = 7rfcG(k) 



7A: + /3' 



27rfc2 



(23) 



7fc + /3 

The pseudo-energy and pseudo-enstrophy are 
quadratic magnitudes and therefore Eqs. (|22l) and 
(|23p must be positive. The relation 7fc > — /3 must apply 
to every value of k, and as the case fc = 1 is the more 
restrictive, we simply need the system to satisfy the 
condition 



(24) 



This condition is enough for integral (|18p to converge. 
Furthermore, asking the total pseudo-energy and pseudo- 
enstrophy to be positive, it is obtained that 7 > and 
that the pseudo-enstrophy thermalizes. 

A truncated inviscid SQG system is expected to reach 
at large times the absolute equilibrium solutions ([^ and 
The values of 7 and /3 are uniquely determined by 
the total amount of pseudo-energy and pseudo-enstrophy 
contained in the initial conditions, and can be calculated 
solving the system of equations 

i?(t = 0)=V-— — , (25) 



^7fc + /3' 



G{t = 0) = ^ 



27rfc2 
7fc + /3' 



(26) 



Note these equations can be easily converted to two 
polynomial equations with two unknowns (7 and (3). 
The number of terms on the r.h.s. of Eqs. (P5|) and 
(|26p depend on the maximum wavenumber preserved af- 
ter the truncation (and therefore, in numerical simula- 
tions, on the linear resolution N). As E{t = 0) and 
G{t = 0) are known, the system can be easily solved 
numerically. For the resolutions used in the next sec- 
tion, the values obtained for (7, /3) are: (346.23, —342.60) 
for N = 32, (1344.79,-1341.39) for N = 64, 
(5044.05,-5040.75) for N = 128, (19504.20,-19500.90) 
for N = 256, (77579.35, -77576.13) for N = 512, and 
(311246.27, -311243.06) for N ^ 1024. 



III. NUMERICAL RESULTS 

In this section we present numerical simulations of the 
SQG equations. We first compare the inviscid spectra 
at late times stemming from the simulations with the 
results derived above. Then, we study the transition from 
the initial condition towards the equilibrium spectrum 
in simulations at larger spatial resolution, and compare 
the behavior of the non-thcrmalized fraction of pseudo- 
energy and pseudo-enstrophy in the inviscid runs with 
the pseudo-energy an pseudo-enstrophy in viscous runs. 
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FIG. 1: Isotropic pseudo-enstrophy spectra for runs with 
= 32 (diamonds), 64 (squares), and 128 (triangles) at 
t = 20, 1000, and 24000 respectively. The Gibbs ensemble 
prediction (|23[) for each run is shown in daslied lines. The 
values of (7, /3) obtained from solving Eqs. (|25p and (|25p for 
the initial pseudo-energy and pseudo-enstrophy of the runs 
are (346.23,-342.60) for iV = 32, (1344.79,-1341.39) for 
iV = 64, and (5044.05, -5040.75) for N = 128. Theoretical 
and numerical results are in good agreement for all wavenum- 
bers. 




FIG. 2: Isotropic pseudo-enstrophy spectra att = 200 for runs 
with — 256 (diamonds), 512 (triangles), and 1024 (crosses). 
The Gibbs ensemble prediction (|23|) is shown in dashed lines. 
The values of (7, (3) obtained from solving Eqs. (|25p and (|25p 
for the initial pseudo-energy and pseudo-enstrophy of the runs 
are (19504.20, -19500.90) for iV = 256, (77579.35^ -77576.13) 
for N = 512, and (311246.27, -311243.06) for iV = 1024. As 
the resolution increases, it takes longer times for the equilib- 
rium at large scales to be reached, while thermalization at 
large wave numbers is achieved fast. 



A. Inviscid truncated runs 

To solve numerically Eq. (|4]) in a 2D periodic domain of 
length 27r, we use a parallel pseudospectral code [47l.l48j. 
Second order Runge-Kutta is used to evolve in time, 
and non- linear terms are computed with the 2/3-rule for 
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FIG. 3: Time for 50% of the pseudo-enstrophy to thermal- 
ize in the inviscid runs as a function of the resohition A^. 



dealiasing, so that the Fourier space is truncated at the 
maximum wave number kmax = N/3, with N the hn- 
ear resolution. Under these conditions, pseudospectral 
methods are known to be equivalent to Galcrkin spectral 
methods (49| , and all quadratic invariants of the original 
equations are conserved in the truncated Fourier space. 
As an example, in simulations of the inviscid SQG equa- 
tions with resolution of 512^ grid points and a time step 
of At = 2 X 10^**, the pseudo-energy was conserved up 
to the fourth digit after 500 turnover times. 

Spatial resolutions rang ed from 32^ to 1024^ grid 
points. The initial condition for all runs is a super- 
position of Fourier modes with random phases, with a 
pseudo-energy spectrum E{k) ~ k~'^ for 1 < fc < 4, 
and zero otherwise. All simulations behave in a similar 
fashion, showing a progressive thermalization of modes 
with higher wavenumber and reaching equilibrium for 
long times; however, simulations at larger resolution take 
longer time to reach equilibrium at all scales. Once the 
system reaches equilibrium, its spectra should be compat- 
ible with solutions and . To compare simulation 
results with the theoretical predictions, 7 and /3 were de- 
termined solving the set of Eqs. (^5)) and (^5]) for each 
spatial resolution separately, as explained in the previous 
section. 

Numerical results for the isotropic pseudo-enstrophy 
spectrum once the equilibrium was reached are shown 
in Fig. [1] with different symbols for the numerical data, 
and with dashed lines for the theoretical predictions. Nu- 
merical spectra are in good agreement with the theory 
at all wavenumbers, and the spectrum shows a peak at 
k = 1 (associated with the condensation of pseudo-energy 
at the gravest modes), and a ^ fc scaling (associated 
with the thermalization of pseudo-enstrophy) for larger 
wavenumbers. However, note that as the resolution is 
increased, the time to reach the equilibrium increases. 

Results for higher resolutions (from 256'^ to 1024^ grid 
points) at t = 200 together with their theoretical pre- 
dictions arc shown separately in Fig. [51 For higher reso- 
lutions, the thermalization process is slower and conver- 
gence of lower wavenumber modes to the statistical equi- 
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FIG. 4: Evohition of the pseudo-energy spectra for the in- 
viscid run with A'^ = 1024 from t = 1 to t = 205 with time 
increments At = 4. A k~^^^ power law followed by an equi- 
librium subrange is observed for early times. As time evolves, 
more modes approach the equihbrium spectrum. The Gibbs 
ensemble prediction is indicated as a reference by the hori- 
zontal line. The arrow indicates the direction in which the 
spectrum evolves in time. 
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FIG. 5: Evolution of the pseudo-enstrophy spectra for the in- 
viscid run with A*' = 1024 from t = 1 to t = 205 with time 
increments At = 4. A k~^^^ power law followed by a ther- 
malized subrange is observed for early times. As time evolves, 
more modes achieve thermalization. The Gibbs ensemble pre- 
diction is indicated as a reference by the ~ k line. The arrow 
indicates the direction in which the spectrum evolves in time. 



librium solution takes progressively longer times. Never- 
theless, numerical and theoretical results agree well for 
intermediate to high wavenumbers. 

The slow down of the thermalization as resolution is 
increased can also be seen in Fig. [3l that shows the time 
t* for the system to have 50% of its pseudo-enstrophy 
thermalized. Note this thermalization takes place pre- 
dominantly at large wavenumbers, and to reach the equi- 
librium at lower wavenumbers takes substantially longer 
times. 

Simulations at lower resolution provide as a result a 
faster way to verify the validity of the Gibbs ensemble 
prediction. However, to study the viscous-like transient 
that develops as the system evolves towards equilibrium. 
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FIG. 6: Time evolution of the thermalized pseudo-energy 
(thin solid) and pseudo-enstrophy (thin dashed line) normal- 
ized by the total pseudo-energy and pseudo-enstrophy for the 
1024'^ run. Both magnitudes grow monotonically towards 
their equilibrium asymptotic value. Non-thermalized pseudo- 
energy (thick solid) and pseudo-enstrophy (thick dashed) 
normalized respectively by total pseudo-energy and pseudo- 
enstrophy are also shown. 



FIG. 7: Time evolution of kth for the inviscid run with 
— 1024. As time evolves, more modes reach thermaliza- 
tion, but kth does not converge to zero as condensation of 
pseudo-energy takes place at the gravest modes. Inset: time 
evolution of Lth ~ 2n/kth for the inviscid runs with A'' — 1024 
(solid) and with A'^ — 512 (dashed) with a ~ t^^* slope shown 
as a reference. 



the runs with larger resolutions will allow us better iden- 
tification of scaling laws and comparison with viscous 
runs. We therefore focus in the following on the 1024^ 
run. Figures U] and [S] show the time evolution of the 
pseudo-energy and pseudo-enstrophy spectra in this run. 
At early times both spectra develop a viscous-like in- 
ertial range, with slopes compatible with Kolmogorov- 
Batchelor-Kraichnan phenomenology [1^. A power law 
^ can be identified in the pseudo-enstrophy, and 

^ k^^/^ in the pseudo-energy spectrum. As time evolves, 
the pseudo-enstrophy shows a progressive thermalization 
starting from the largest wavenumbers, and the viscous- 
like inertial range becomes narrower as the G{k) ~ k 
thermalized spectrum broadens. We will take the mini- 
mum of the pseudo-enstrophy spectrum, at fc = kth, as 
the delimiting wavenumber between these two subranges. 
The flat pseudo-energy spectrum for k > kth, and its 
peak at fc = 1, evidences a condensation of this quantity 
at low k rather than a thermalization at high wavenum- 
bers. 

In the following we define the thermalized pseudo- 
energy and pseudo-enstrophy respectively as the sum of 
the pseudo-energy and pseudo-enstrophy contained in all 
modes with fc > kth, 

Eth = 2 ^(^)' = £ ^(^)- (27) 

k—kih k—kth 

In a similar way, we define the non-thermalized pseudo- 
energy and pseudo-enstrophy as the difference between 
the total amount of these quantities (which is constant 
in time) and their thermalized values, 

Enth — E — Eth, Gnth ~ G — Gth- (28) 



The time evolution of Eth, Enth, Gth, and Gnth nor- 
malized respectively by E and G is shown in Fig. |6l 
Eth and Gth grow monotonically at early times, converg- 
ing to an almost constant value for long times. This 
is consistent with Figs. U and [SJ where more and more 
modes approach the statistical equilibrium solution as 
time evolves. This is also illustrated by the evolution of 
kth in Fig. [T] As most of the pseudo-energy condenses at 
fc = 1, Eth remains a small fraction of the total pseudo- 
energy while Enth (which in this case represents the con- 
densed pseudo-energy) stays almost constant. On the 
other hand, Gth grows to a larger fraction of the total 
pseudo-enstrophy as the ~ fc thermalized spectrum for 
G at large wavenumbers (see Fig. [5|) holds a significant 
fraction of the pseudo-enstrophy in the system. 

As shown in Ref. [3| for the truncated Euler equations, 
although the truncated system is inviscid, the transfer of 
pseudo-enstrophy towards thermalized modes before the 
equilibrium is reached (resulting in the growth of Eth and 
Gth in time) can be interpreted as a viscous-like transient 
in which the thermalized modes give rise to an effective 
viscosity acting on the non-thermalized range. This ef- 
fective viscosity is responsible for the development of tur- 
bulent inertial subranges. In the next subsection we ana- 
lyze in more detail this transient, comparing the inviscid 
system with viscous solutions. 



B. Inviscid vs. dissipative systems 

We focus now on the period of time when Eth and Gth 
are time dependent, and the flux of pseudo-enstrophy 
towards thermalized modes can be interpreted as an out- 
of-cquilibrium turbulent solution for k < kth, with ef- 
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FIG. 8: Pseudo-energy spectra at different times for 1024'' in- 
viscid (solid) and viscous (dasfied) runs. Botli runs develop a 
power law compatible with ~ k~^^^ (indicated as a reference 
by the straight line). Later on, the spectra differ at inter- 
mediate and large wavenumbers, with the gravest modes still 
showing similar amplitudes. The arrows indicate the values 
of kth at each time. 



FIG. 9: Pseudo-enstrophy spectra at different times for the 
same 1024^ inviscid (solid) and viscous (dashed) runs shown 
in Fig. [SI A ~ k~^^'^ slope is shown as a reference. The arrows 
indicate the values of kth at each time. 
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fective viscosity associated with the thermahzed modes 
with k> kth- We show that during this period, the ideal 
truncated model can give valuable information on the be- 
havior of a similar but dissipative system (by similar, we 
mean a viscous system subject to the same initial condi- 
tions). With this aim, wc compare ideal and dissipative 
spectra and decays for this particular period of time. 

To consider dissipative SQG flows, wc must solve 
Eq. (|4]) with the addition of a dissipative term. 



^ = ( 5y-0 - dxi^ dyq 



:/|k|V- 



(29) 



We solve this equation in the same way as in the ideal 
case, with a pseudospcctral method with the 2/3-rule for 
dealiasing, and with second order Runge-Kutta to evolve 
in time. 

Equation (p9|) in terms of the scalar q in real space 
reads 



dt 



-u-\7q + vV^q 



(30) 



Multiplying the equation by q the following conservation 
law follows 



2~dt 



1 



V • (q^u) + lyqV^q. 



(31) 



Integrating over the entire domain, the first term on the 
r.h.s. cancels and Eq. ([31]) in Fourier space results 



— = — / G{k)dk = -2v / k^G{k)dk = -o- (32) 
dt dt Jq Jq 

with a the pscudo-cnstrophy dissipation rate. Introduc- 
ing a dissipation wavenumber /c,, such that all dissipa- 
tion is concentrated in Fourier space between k = I and 
k = k^i, in the turbulent steady state we can approximate 



cr « 2i^ 



k^G{k)dk. 



(33) 



From Ref. [Hi, G{k) ~ a^/^fc-^/^, and replacing in 
Eq. (I2S1) wc finally get 



(34) 



This relation is important to fix the resolution in vis- 
cous simulations so that all relevant scales in the flow 
are well resolved. In practice, we want the dissipation 
wavenumber to be smaller than the maximum resolved 
wavenumber kmax, or in other words, we ask for the con- 
dition 



E 



k^G{k) 



1/4 



< 



N 



to be fulfilled at all times. 






100.0 



time 



FIG. 10: Decay of pseudo-enstrophy G(t) in viscous runs 
(solid line) and time evolution of Gnth(i) = G — Gth(t) in 
ideal truncated runs (dashed lines) at increasing resolution: 
(a) iV = 128, u 10~^ (b) iV = 256, i/ = 5 x 10"*, (c) 

TV = 512, = 2.5 X 10""*, and (d) N = 1024, v = x 10"^ 
The viscous G(t) approaches Gnthif) as the Reynolds number 
is increased. A ~ t"^/* decay is indicated as a reference. The 
(35) dissipation wavenumbers for each viscous run are: kr, = 36 



150 for N = 512, 



for iV = 128, fc„ = 84 for = 256, h 
and krf = 318 for N = 1024, while the maximum resolved 
wavenumber is kmax ~ N/3. 
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FIG. 11: Decay of pseudo-enstrophy G{t) for viscous runs 
iV = 1024, N = 512, iV = 256, A*' = 128 from top to bottom 
in semi-log scale. Note the exponential decay after t ~ 80. 



In Fig. IS] we show the pseudo-energy spectra at dif- 
ferent times for the 1024^ ideal truncated run, and for 
a 1024^ well resolved viscous run with i' = 4 x 10^^. 
At early times, the pseudo-enstrophy contained at large 
scales is transferred to smaller scales, and is eventually 
affected by thermalization or by viscous dissipation in 
the ideal or viscous case respectively. These effects are 
different, as thermalization is a time dependent problem 
which results in a time dependent effective viscosity (see, 
e.g., O [Hi), while viscosity in the dissipative run is 
fixed. However, the large scales of both systems are sim- 
ilar at low wavenumbers, with both systems developing 
a ^ k~^/'^ power law in the pseudo-energy in accordance 
with the theoretical results (although this subrange is 
wider at early times in the inviscid truncated run). This 
power law is lost at later times in the inviscid run, al- 
though the amplitude of the gravest modes still shows 
agreement with the viscous run. A similar behavior is 
observed in the time evolution of the pseudo-enstrophy 
spectrum (see Fig. [S]). 

The similitude in the evolution of the large scale spec- 
tra in viscous and inviscid runs makes it plausible to 
compare the decay of quadratic quantities in the viscous 
simulations, say G(i), with the evolution of the associ- 
ated non-thermalized quantity Gnth{i) in inviscid trun- 
cated runs. In freely decaying turbulent flows, quadratic 
quantities often develop a self-similar decay law in time, 
once turbulence develops and dissipation sets in. The 
power law obeyed during this decay (in this example, say, 
G{t) t~^) often requires large resolution as viscous ef- 
fects tend to affect the decay making determination of 
the decay rate from the simulations difficult. 

In Fig. [ini we compare the decay of the pseudo- 
enstrophy in several dissipative runs with the evolution 
of Gnth{t) in inviscid runs at different grid sizes. In all 
viscous cases, the pseudo-enstrophy is quasi-conserved at 
early times, and then a self-similar decay develops fol- 
lowing a power law compatible (at the largest Reynolds 



FIG. 12: Evolution of G„th for different grid sizes: = 128 
(dashed), 256 (dash-dotted), 512 (dash-triple-dotted), and 
1024 (dotted). There is small delay in the time when the 
decay begins as the grid size (or equivalently, kmax) is in- 
creased, shown in more detail in the inset. Except for this 
delay, all runs show the same decay, with a ~ t~^^'^ decay 
shown as a reference. 

numbers and spatial resolution) with G{t) t''^/^. It 
is worth noting that such a decay is also compatible 
with strict mathematical bounds for the decay of pseudo- 
enstrophy derived in [H, . The decay is followed by a 
saturation and a viscous exponential decay of the pseudo- 
enstrophy that remains in the system (see Fig. [TT|| . 

The non-thcrmalizcd Gnth{t) in all the inviscid runs 
shows the same evolution, which except for a small delay 
in the time when the quasi-conservation is broken (see 
Fig. [T^ and inset), is independent of the spatial resolu- 
tion considered. Remarkably, the evolution of G{t) seems 
to approach, as the Reynolds number is increased, the be- 
havior of Gnthit), to the point that the decay t~^^^ can 
be identified in Gnthit) at much lower resolution (e.g., 
grid sizes of 128^ or 256^ grid points) that in the viscous 
run. 

In this light, the decay of viscous SQG can be under- 
stood as follows. At early times, G(t) remains approxi- 
mately constant as the inertial range spectrum develops. 
As the smallest scales have not been excited yet, dissipa- 
tive effects are negligible. Once turbulence develops, the 
decay law G{t) ~ t~^^^ is observed. The decay rate in 
this stage is dominated by the turbulent cascade of the 
pseudo-enstrophy towards smaller scales, and the time 
evolution is given by the balance equation dG/dt = —a 
where the flux of pseudo-enstrophy a is controlled by the 
non-linear term in the equation, and independent of the 
value of the viscosity v as long as the Reynolds number 
is large enough. As a result, the truncated inviscid run, 
although it has a different effective viscosity, shows the 
same decay. Later, all pseudo-enstrophy that remains in 
the system is associated with the condensation of pseudo- 
energy at fc = 1, and this remaining pseudo-energy (as 
well as the pseudo-enstrophy) can only decay exponen- 
tially in the viscous run fFig. Ill|). while it remains con- 
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FIG. 13: Comparison between E„th in an TV = 1024 trun- 
cated inviscid simulation (dashed line) and E{t) for viscous 

256 {u = 5x10-''), 



simulations with = 128 (i^ = 2 x 10 

-4 



), and 1024 (i/ = 9 x lO"'') (bottom to 



512 (i/ = 2.5 X 10" 
top). 



stant in the inviscid run (Fig. HZ 

Figure lT^ compares the non-thermalized fraction of the 
pseudo-energy Enth{t) = E - Eth{t) in the N = 1024 
simulation with the decay of E{t) in the viscous simula- 
tions with N = 128, 256, 512, and 1024. In the inviscid 
case, as most of the pseudo-energy condenses, there is no 
significant change in E^th with time. The viscous runs 
approach this behavior as the viscosity is decreased, in 
agreement with Batchelor-Kraichnan-Leith phenomenol- 
ogy for systems with condensation of one invariant at 
large scales (see, e.g., similar arguments for the decay of 
the enstrophy while the energ y r emains approximately 
constant in 2D Navier-Stokes [50] )• It is worth noting 
that a slow down in the decay of the pseudo-energ y a s 
the viscosity was decreased was already observed in [43| ■ 

So far, random initial conditions were used in all runs. 
To analyze sensitivity of the decay to other initial condi- 
tions, we briefly discuss results for three initial conditions 
often used in studies of singularities in the QG equations 



I) q{t = 0) = sin(.T) sin(y) -|- cos(y). 



(36) 



II) q{t = 0) = -(cos(2x) cos{y) + sin(x) sin(2/)), (37) 

III) q{t = 0) = cos(2x) cos(j/) -I- sin(x) sin(j/) -|- 

-Kcos(2a;)sin(3?/). (38) 

Initial conditions II and III lead to the same decay of 
the pseudo-enstrophy ^ t^^/^ as observed in the previous 
simulations (see Fig. fH)) . However, initial condition I re- 
sults in a much slower decay (not shown), and up to the 
time we integrated the equations we were unable to iden- 
tify any power law. Further analysis showed that the flow 
resulting from this initial condition, which initially acti- 
vates only a few modes in Fourier space, develops a much 



FIG. 14: Evolution of the pseudo-enstrophy for simulations 
with initial conditions II and III, = 1024 and u = Sjx 10"^. 
Both runs show a decay compatible with G{t) ~ t~^^'^ . 



steeper pseudo-energy and pseudo-enstrophy spectra. As 
discussed above, a fully developed turbulent spectrum is 
needed to observe the self-similar decay in the enstrophy. 

As previously mentioned, the G{t) ^ t~^^^ decay ob- 
served in the simulations is consistent with the upper 
limit of strict mathematical bounds derived for the de- 
cay of pseudo-enstrophy (see, e.g., [H, EO])- From the 
physical point of view, it is interesting that such a decay 
law can be related with the decay followed by other sys- 
tems that develop condensates in the inviscid truncated 
case, and large-scale coherent structures in the viscous 
case (see, e.g., [1^). In 2D Navier-Stokes turbulence, 
the separation of scales between energy containing scales 
and enstrophy containing scales (as the former is con- 
centrated at large scales while the latter is concentrated 
at small scales), together with assumptions of weak cou- 
pling between these scales [H^ [s^l , led to the idea that 
the decay of 2D Navier-Stokes turbulence can be under- 
stood as the superposition of a coherent (condcnsated) 
state, and a random phase, separated by a "fuzzy scale" 
following ideas developed for Bose-Einstein condensates 
[ssf . In SQG, from the results in SecHH we can consider 
a pseudo-energy containing scale associated with the con- 
densation of pseudo-energy at low wavenumbers, which 
contains almost all of the system pseudo-energy, and a 
pseudo-enstrophy containing scale concentrating most of 
the pseudo-enstrophy. In the inviscid truncated case, 
we can easily identify the fuzzy scale delimiting these 
two phases with the thermalization length Lth = 'i.'n /ktu- 
From Eq. ([5^ . dimensional analysis leads to 



dG 
'dt 



0^2 
L 



(39) 



where L is a characteristic scale. Associating L = Lth, 
as most of the pseudo-enstrophy is contained below that 
scale, Eq. (jJH) is only consistent with the observed decay 
G{t) - if Lthit) - f^^. The inset in Fig. [TU] shows 

the time evolution of Lth in the truncated inviscid runs 
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with N = 1024 and A'' = 512 with random initial condi- 
tions, appears a reasonable fit; the other simulations 
are also consistent with this power law. 

IV. CONCLUSIONS 

We derived statistical equilibrium solutions of the trun- 
cated inviscid SQG equations, and verified their validity 
at late times in numerical simulations of the truncated 
SQG equations. Numerical spectra arc in agreement with 
the theory at all wavenumbers, although as resolution is 
increased it takes longer times for the system to reach 
equilibrium. 

We also studied the evolution towards equilibrium 
of the ideal runs, finding that both pseudo-energy and 
pseudo-enstrophy spectra develop a viscous-like inertial 
range with slopes ^ k~^/^ and ~ k~^^^ respectively, and 
compatible with Kolmogorov-Batchelor-Kraichnan phe- 
nomenology. Note that in previous studies of the invis- 
cid SQG equations, a steeper power law was observed 
[4^ . but the initial conditions I that were used result 
in steeper spectra for the pseudo-energy and pseudo- 
enstrophy. As time evolves, the pseudo-enstrophy in our 
runs shows a gradual thermalization of higher wavenum- 
ber modes with G(k) ~ k, narrowing the viscous-like 
inertial range. The spectrum for the pseudo-energy 
presents a flat scaling for high wavenumbers and a peak 
at fc = 1, showing condensation of pseudo-energy at small 
wavenumbers instead of thermalization. 

Defining Gth and Eth respectively as the pseudo- 
enstrophy and pseudo-energy contained in the thermal- 
ized modes, we studied the period of time during which 
these quantities are time-dependent. Through this pe- 
riod, the transfer of pseudo-enstrophy towards thermal- 
ized modes results in a viscous-like cascade with effective 
viscosity associated with the thermalized modes. This al- 
lowed us to compare spectra and decay in ideal runs with 
runs subjected to identical initial conditions but with the 
addition of dissipation. The large scales of both systems 
arc similar and both systems develop a ~ k~^/^ power 
law for the pseudo-energy in accordance with the theo- 
retical results. At later times, the power law is lost due 
to thermalization in the inviscid runs and dissipation in 
the viscous runs. 

We compared the free decay of G{t) with the time evo- 
lution of the non-thermalized Gnth{t) in inviscid trun- 



cated runs for different grid sizes. Following an initial 
period of time in which the pseudo-enstrophy is quasi- 
conserved, viscous cases develop a power-law decay com- 
patible with G{t) ~ for the largest Reynolds num- 
bers and spatial resolutions considered. Later in time, 
the decay becomes shallower and an exponential viscous 
decay of the remaining pseudo-enstrophy follows. The 
power-law decay is in agreement with strict mathemati- 
cal bounds for the decay of the pseudo-enstropy [H, , 
and seems to be robust for initial conditions in which 
most of the pseudo-energy is held at the lowest modes. 
From the point of view of the statistical equilibrium so- 
lutions, the decay law can be understood as the result 
of two decaying phases weakly interacting and separated 
by a fuzzy scale: a coherent (condensated) state that 
contains most of the pseudo-energy, and a random phase 
that contains most of the pseudo-enstrophy. 

With the exception of a small delay in the time when 
the quasi-conservation is broken, all the inviscid runs 
show the same evolution for the non-thermalized Gnth (t) 
independently of the spatial resolution considered. The 
evolution of G{t) approaches that of Gnth as the Reynolds 
number is increased. Remarkably, inviscid systems as 
small as 128^ already give information about the free de- 
cay law expected for high Reynolds viscous systems. 

This good agreement between the evolution of the non- 
thermalized components of the pseudo-enstrophy and the 
pseudo-energy in truncated inviscid runs at low reso- 
lution, and of the decay of pseudo-energy and pseudo- 
enstrophy in viscous runs at large resolution, further in- 
dicate a new application of inviscid simulations: that of 
estimating the decay law followed by the viscous system 
at very large Reynolds number. In this light, the decay 
of pseudo-enstrophy in SGQ turbulence in the limit of 
very large Reynolds number would be compatible with 
G{t) ^ t~^/^ when all initial perturbations are at the 
gravest modes in the system. It is unclear for the mo- 
ment whether the evolution of non-thermalized quanti- 
ties in other inviscid truncated system can be also used 
as proxies of decay. This problem is left for future work. 
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